```{r}

rm(list = ls())
library(pacman)
p_load("fixest", "haven", "conleyreg", "readxl", "dplyr", "lmtest", "ggplot2", "ggthemes", "jtools", "tidyverse", "scales")

Panel_long <- read_xlsx("C:/Users/tanay/Dropbox/ADS/Panel_28Jan2024.xlsx")


```


#Main Regressions
```{r}

Regressions_main <- list(

##TWFE  
    ##Human deaths
  
feols(HDI_scaled ~ ADS2 + ADS_IN_0_1_2  + neigh + dist | code + year, data = Panel_long, conley(25)),

   ##Elephant deaths

feols(EDI_scaled ~ ADS2 + ADS_IN_0_1_2  + neigh +  HD_dummy_0_1 + dist | code + year, data = Panel_long, conley(25)),

##Year FE  
    ##Human deaths
  
feols(HDI_scaled ~ ADS_ever + ADS2 + ADS_IN_0_1_2 + ADS_nbd_ever + neigh + dist| year, data = Panel_long, conley(25)),

   ##Elephant deaths

feols(EDI_scaled ~ ADS_ever + ADS2 + ADS_IN_0_1_2  + ADS_nbd_ever + neigh +  HD_dummy_0_1 + dist | year, data = Panel_long, conley(25)),

##Poisson
  ##Human deaths

fepois(HDI ~ ADS2 + ADS_IN_0_1_2  + neigh + dist | code + year, data = Panel_long, cluster = "code"),

  ##Elephant deaths
fepois(EDI ~ ADS2 + ADS_IN_0_1_2  + neigh + HD_dummy_0_1 + dist | code + year, data = Panel_long, cluster = "code"))


## Getting coefficients and robust variances of Poisson regressions

Objects <- list(
  HDI = tibble(Varname =  rownames(coeftest(Regressions_main[[5]])), 
                  Coeff = coeftest(Regressions_main[[5]])[, 1],
                  Var = diag(vcov_cluster(Regressions_main[[5]], "code"))),
  
  EDI = tibble(Varname =  rownames(coeftest(Regressions_main[[6]])), 
                  Coeff = coeftest(Regressions_main[[6]])[, 1],
                  Var = diag(vcov_cluster(Regressions_main[[6]], "code"))))

##Calculating IRR = e^b and its se = sqrt(e^b*var*e^b) = IRR*se

Objects <- Map(function(x) {
    mutate(x, IRR = exp(Coeff),
           )}, Objects)

Objects <- Map(function(x) {
    mutate(x, IRR_se = IRR*sqrt(Var),
           )}, Objects)

Objects <- Map(function(x) {
    mutate(x, IRR_z = IRR/IRR_se,
           )}, Objects)


Regressions_main[[5]]$coefficients <- (as.vector(Objects$HDI[, 4]))[[1]] 
Regressions_main[[5]]$se <- (as.vector(Objects$HDI[, 5]))[[1]] 


Regressions_main[[6]]$coefficients <- (as.vector(Objects$EDI[, 4]))[[1]]
Regressions_main[[6]]$se <- (as.vector(Objects$EDI[, 5]))[[1]]

##Constructing VCOV matrices for regressions

Vcov <- list(
  vcov_conley(Regressions_main[[1]], cutoff = 25),
  vcov_conley(Regressions_main[[2]], cutoff = 25),
  vcov_conley(Regressions_main[[3]], cutoff = 25),
  vcov_conley(Regressions_main[[4]], cutoff = 25),
  (diag(exp(c(Objects$HDI$Coeff))))%*%vcov_cluster(Regressions_main[[5]], cluster = "code")%*%(diag(exp(c(Objects$HDI$Coeff)))),
  (diag(exp(c(Objects$EDI$Coeff))))%*%vcov_cluster(Regressions_main[[6]], cluster = "code")%*%(diag(exp(c(Objects$EDI$Coeff))))) 

outcome_vars <- c("HDI", "EDI", "HDI_scaled", "EDI_scaled", "EDI_accident", "EDI_accident_scaled", "EDI_other", "EDI_other_scaled")

Baseline_avg_012 <- filter(Panel_long, ADS_ever == 1 & ADS_IN_0_1_2 == 0 & ADS2 == 0) %>% dplyr::select(all_of(outcome_vars)) %>% lapply(function(x) mean(x))

Baseline_avg_01 <- filter(Panel_long, ADS_ever == 1 & ADS_IN_0_1 == 0 & ADS2 == 0) %>% dplyr::select(all_of(outcome_vars)) %>% lapply(function(x) mean(x))

Baseline_avg_0 <- filter(Panel_long, ADS_ever == 1 & ADS_IN_0 == 0 & ADS2 == 0) %>% dplyr::select(all_of(outcome_vars)) %>% lapply(function(x) mean(x))


##Table: Main Regressions
etable(list(Regressions_main[[1]], Regressions_main[[2]]), title = "Effects of Anti-Depredation Squads on mortality due to conflict", vcov = list(Vcov[[1]], Vcov[[2]]), fitstat = c("n", "r2", "aic", "bic"), coefstat = "se", style.tex = style.tex(main = "base"), replace = T, drop.section = c("fixef"), signif.code =  c("***"=0.01, "**"=0.05, "*"=0.10), extralines = list("Baseline" = c(Baseline_avg_012$HDI_scaled, Baseline_avg_012$EDI_scaled), "Fixed effects" = c("Two-way", "Two-way"), "Standard Errors" = c("Conley", "Conley")), dict = c("Pop" = "Population", "ADS2" = "ADS", "dist" = "dist_eleph_habitat", "neigh" = "ADS_neighbor", "HDI_scaled" = "Human deaths", "EDI_scaled" = "Elephant deaths", "HD_dummy_0_1" = "HD_0_1"), file = "C:/Users/tanay/downloads/Main_Regressisons_LM_TWFE.tex")


etable(list(Regressions_main[[3]], Regressions_main[[4]], Regressions_main[[5]], Regressions_main[[6]]), title = "Effects of Anti-Depredation Squads on mortality due to conflict", vcov = list(Vcov[[3]], Vcov[[4]], Vcov[[5]], Vcov[[6]]), fitstat = c("n", "r2"), coefstat = "se", style.tex = style.tex(main = "base"), replace = T, drop.section = c("fixef"), signif.code =  c("***"=0.01, "**"=0.05, "*"=0.10), extralines = list("Baseline" = c(Baseline_avg_012$HDI_scaled, Baseline_avg_012$EDI_scaled, Baseline_avg_012$HDI, Baseline_avg_012$EDI), "Fixed effects" = c("Year", "Year", "Two-way", "Two-way"), "Standard Errors" = c("Conley", "Conley", "Clustered", "Clustered")), dict = c("Pop" = "Population", "ADS2" = "ADS", "dist" = "dist_eleph_habitat", "neigh" = "ADS_neighbor", "HDI_scaled" = "Human deaths", "EDI_scaled" = "Elephant deaths", "HDI" = "Human deaths", "EDI" = "Elephant deaths", "HD_dummy_0_1" = "HD_0_1"), label = "tbl:main_YFE_Poi", file = "C:/Users/tanay/downloads/Main_Regressisons_Poisson_YFE.tex")

```


##EDI regressions

```{r}

Regressions2 <- list(

## accident - OLS

feols(EDI_accident_scaled ~ ADS2 + ADS_IN_0_1_2  + neigh +  HD_dummy_0_1 + dist | code + year, data = Panel_long, conley(25)),

## other - OLS

feols(EDI_other_scaled ~ ADS2 + ADS_IN_0_1_2  + neigh +  HD_dummy_0_1 + dist | code + year, data = Panel_long, conley(25)),

## accident - Poisson

fepois(EDI_accident ~ ADS2 + ADS_IN_0_1_2  + neigh + HD_dummy_0_1 + dist | code + year, data = Panel_long, cluster = "code"),

## other - Poisson

fepois(EDI_other ~ ADS2 + ADS_IN_0_1_2  + neigh + HD_dummy_0_1 + dist | code + year, data = Panel_long, cluster = "code"))


## Getting coefficients and robust variances of Poisson regressions

Objects <- list(
  EDI_accident = tibble(Varname =  rownames(coeftest(Regressions2[[3]])), 
                  Coeff = coeftest(Regressions2[[3]])[, 1],
                  Var = diag(vcov_cluster(Regressions2[[3]], "code"))),
 
  EDI_other = tibble(Varname =  rownames(coeftest(Regressions2[[4]])), 
                  Coeff = coeftest(Regressions2[[4]])[, 1],
                  Var = diag(vcov_cluster(Regressions2[[4]], "code"))))


##Calculating IRR = e^b and its se = sqrt(e^b*var*e^b) = IRR*se

Objects <- Map(function(x) {
    mutate(x, IRR = exp(Coeff),
           )}, Objects)

Objects <- Map(function(x) {
    mutate(x, IRR_se = IRR*sqrt(Var),
           )}, Objects)

Objects <- Map(function(x) {
    mutate(x, IRR_z = IRR/IRR_se,
           )}, Objects)

Regressions2[[3]]$coefficients <- (as.vector(Objects$EDI_accident[, 4]))[[1]]
Regressions2[[3]]$se <- (as.vector(Objects$EDI_accident[, 5]))[[1]]

Regressions2[[4]]$coefficients <- (as.vector(Objects$EDI_other[, 4]))[[1]]
Regressions2[[4]]$se <- (as.vector(Objects$EDI_other[, 5]))[[1]]


##Constructing VCOV matrices for regressions

Vcov <- list(
  vcov_conley(Regressions2[[1]], cutoff = 25),
  
  vcov_conley(Regressions2[[3]], cutoff = 25),
  
  (diag(exp(c(Objects$EDI_accident$Coeff))))%*%vcov_cluster(Regressions2[[3]], cluster = "code")%*%(diag(exp(c(Objects$EDI_accident$Coeff)))),
    
  
  (diag(exp(c(Objects$EDI_other$Coeff))))%*%vcov_cluster(Regressions2[[3]], cluster = "code")%*%(diag(exp(c(Objects$EDI_other$Coeff))))) 


##Table: EDI
etable(Regressions2, title = "Effects of Anti-Depredation Squads on elephant deaths (by type)", vcov = Vcov, fitstat = c("n", "aic", "bic"), style.tex = style.tex(main = "base"), replace = T, signif.code = c("***"=0.01, "**"=0.05, "*"=0.10), drop.section = c("fixef"), extralines = list("Baseline" = c(Baseline_avg_012$EDI_accident_scaled, Baseline_avg_012$EDI_other_scaled, Baseline_avg_012$EDI_accident,  Baseline_avg_012$EDI_other), "Fixed effects" = c("Two-way", "Two-way", "Two-way", "Two-way"), "Standard Errors" = c("Conley", "Conley", "Clustered", "Clustered")), dict = c("Pop" = "Population", "ADS2" = "ADS", "dist" = "dist_eleph_habitat", "neigh" = "ADS_neighbor", "HDI_scaled" = "HD", "EDI.wo.Hits_scaled" = "ED w/o hits", "EDI.wo.Hits" = "ED w/o hits", "HD_dummy_0_1" = "HD_0_1", "EDI_accident_scaled" = "EDI (accident)", "EDI_accident" = "EDI (accident)", "EDI_other_scaled" = "EDI (other)", "EDI_other" = "EDI (other)"), file = "C:/Users/tanay/downloads/EDI_bytype.tex")

```

### Figure: Regression summary plot

```{r}

Regs <- list(
  "HD" = tidy(Regressions_main[[1]])[, 1:3], 
  "ED" = tidy(Regressions_main[[2]])[, 1:3],
  "ED (potential accident)" = tidy(Regressions2[[1]])[, 1:3],
  "ED (other)" =  tidy(Regressions2[[2]])[, 1:3]
)



##Figure: linear model 
LM_summary <- plot_summs(Regressions_main[[1]], Regressions_main[[2]], Regressions2[[1]], Regressions2[[2]], colors = c("red", "blue", "blue","blue"), coefs = c("ADS" = "ADS2", "Human death in current or last year" = "HD_dummy_0_1", "ADS introduced within the next two years" = "ADS_IN_0_1_2", "ADS in neighboring village" = "neigh"), line.size = 1, ci_level = 0.95, caption = "Effect of ADS on mortality", model.names = c("Human deaths", "Elephant deaths", "ED: potentially \n accidental", "ED: other"), legend.title = "Outcome") + scale_y_discrete(labels = label_wrap(14))  +  scale_x_continuous(breaks = c(-4, -3, -2, -1, 0, 1, 2, 3, 4)) + theme(axis.text.x=element_text(size=18), axis.text.y=element_text(size=14), legend.text = element_text(size = 14), axis.title.x = element_text(size=14)) 

rm(Regressions2, Regressions_main, LM_summary, Objects, Regs, Vcov)
```



#Robustness: with ADS_IN_0_1 
```{r}

Regressions_r1 <- list(

##TWFE  
    ##Human deaths
  
feols(HDI_scaled ~ ADS2 + ADS_IN_0_1  + neigh + dist | code + year, data = Panel_long, conley(25)),

   ##Elephant deaths

feols(EDI_scaled ~ ADS2 + ADS_IN_0_1  + neigh +  HD_dummy_0_1 + dist | code + year, data = Panel_long, conley(25)),

## ED Regressions 
   ## accident - OLS

feols(EDI_accident_scaled ~ ADS2 + ADS_IN_0_1  + neigh +  HD_dummy_0_1 + dist | code + year, data = Panel_long, conley(25)),

   ## other - OLS

feols(EDI_other_scaled ~ ADS2 + ADS_IN_0_1  + neigh +  HD_dummy_0_1 + dist | code + year, data = Panel_long, conley(25)))


##Constructing VCOV matrices for regressions

Vcov <- list(
  vcov_conley(Regressions_r1[[1]], cutoff = 25),
  vcov_conley(Regressions_r1[[2]], cutoff = 25),
  vcov_conley(Regressions_r1[[3]], cutoff = 25),
  vcov_conley(Regressions_r1[[4]], cutoff = 25)) 


Baseline_avg_01 <- filter(Panel_long, ADS_IN_0_1 == 0, ADS2 == 0) %>% select(outcome_vars) %>% lapply(function(x) mean(x, na.rm = T))

##Table: Main Regressions
etable(list(Regressions_r1[[1]], Regressions_r1[[2]], Regressions_r1[[3]], Regressions_r1[[4]]), title = "Effects of Anti-Depredation Squads on mortality due to conflict", vcov = Vcov, fitstat = c("n", "aic", "bic"), coefstat = "se", style.tex = style.tex(main = "base"), replace = T, drop.section = c("fixef"), signif.code =  c("***"=0.01, "**"=0.05, "*"=0.10), extralines = list("Baseline" = c(Baseline_avg_01$HDI_scaled, Baseline_avg_01$EDI_scaled, Baseline_avg_01$EDI_accident_scaled, Baseline_avg_01$EDI_other_scaled), "Fixed effects" = c("Two-way", "Two-way", "Two-way", "Two-way"), "Standard Errors" = c("Conley", "Conley", "Conley", "Conley")), dict = c("Pop" = "Population", "ADS2" = "ADS", "dist" = "dist_eleph_habitat", "neigh" = "ADS_neighbor", "HDI_scaled" = "Human deaths", "EDI_scaled" = "Elephant deaths", "HD_dummy_0_1" = "HD_0_1"), file = "C:/Users/91941/Desktop/Raghav/CECFEE/Elephant_Paper/Regression Tables/Feb10_24/Regressisons_LM_TWFE_ADS01.tex")


rm(Regressions_r1)
```


#Robustness: w/ ADS_IN_0
```{r}

Regressions_r2 <- list(

##TWFE  
    ##Human deaths
  
feols(HDI_scaled ~ ADS_IN_0 + ADS2  + neigh + dist | code + year, data = Panel_long, conley(25)),

   ##Elephant deaths

feols(EDI_scaled ~ ADS_IN_0 + ADS2 + neigh +  HD_dummy_0_1 + dist | code + year, data = Panel_long, conley(25)),

## ED Regressions

    ## accident - OLS

feols(EDI_accident_scaled ~ ADS_IN_0 + ADS2 + neigh +  HD_dummy_0_1 + dist | code + year, data = Panel_long, conley(25)),

    ## other - OLS

feols(EDI_other_scaled ~ ADS_IN_0 + ADS2 + neigh +  HD_dummy_0_1 + dist | code + year, data = Panel_long, conley(25)))


##Constructing VCOV matrices for regressions

Vcov <- list(
  vcov_conley(Regressions_r2[[1]], cutoff = 25),
  vcov_conley(Regressions_r2[[2]], cutoff = 25),
  vcov_conley(Regressions_r2[[3]], cutoff = 25),
  vcov_conley(Regressions_r2[[4]], cutoff = 25)) 


Baseline_avg_0 <- filter(Panel_long, ADS_IN_0 == 0, ADS2 == 0) %>% select(outcome_vars) %>% lapply(function(x) mean(x, na.rm = T))

##Table: Main Regressions
etable(Regressions_r2, title = "Effects of Anti-Depredation Squads on mortality due to conflict", vcov = Vcov, fitstat = c("n", "aic", "bic"), coefstat = "se", style.tex = style.tex(main = "base"), replace = T, drop.section = c("fixef"), signif.code =  c("***"=0.01, "**"=0.05, "*"=0.10), extralines = list("Baseline" = c(Baseline_avg_0$HDI_scaled, Baseline_avg_0$EDI_scaled, Baseline_avg_0$EDI_accident_scaled, Baseline_avg_0$EDI_other_scaled), "Fixed effects" = c("Two-way", "Two-way", "Two-way", "Two-way"), "Standard Errors" = c("Conley", "Conley", "Conley", "Conley")), dict = c("Pop" = "Population", "ADS2" = "ADS", "dist" = "dist_eleph_habitat", "neigh" = "ADS_neighbor", "HDI_scaled" = "Human deaths", "EDI_scaled" = "Elephant deaths", "HD_dummy_0_1" = "HD_0_1"), file = "C:/Users/91941/Desktop/Raghav/CECFEE/Elephant_Paper/Regression Tables/Feb10_24/Main_Regressisons_LM_TWFE_ADS0.tex")



rm(Regressions2, Regressions_r2)
```


#Robustenss: Additonal Regressors

```{r}

Regressions <- list(

    ##Human deaths
  
feols(HDI_scaled ~ ADS2 + ADS_IN_0_1_2  + neigh +  Pop + dist + prop.Eleph250 + prop.Agri.within  + meanNL  | code + year, data = Panel_long, conley(25)),

   ##Elephant deaths

feols(EDI_scaled ~ ADS2 + ADS_IN_0_1_2  + neigh + HD_dummy_0_1 +  Pop + dist + prop.Eleph250 + prop.Agri.within + meanNL | code + year, data = Panel_long, conley(25)),

##Year fixed effects only
    ##Human deaths
  
feols(HDI_scaled ~ ADS2 + ADS_ever + ADS_IN_0_1_2 + ADS_nbd_ever + neigh +  Pop + dist + prop.Eleph250 + prop.Agri.within  + meanNL  | year, data = Panel_long, conley(25)),

   ##Elephant deaths

feols(EDI_scaled ~ ADS2 + ADS_ever + ADS_IN_0_1_2  + ADS_nbd_ever + neigh + HD_dummy_0_1 +  Pop + dist + prop.Eleph250 + prop.Agri.within + meanNL | year, data = Panel_long, conley(25)),

  
  ##Human deaths

fepois(HDI ~ ADS2 + ADS_IN_0_1_2  + neigh + Pop + dist + prop.Eleph250 + prop.Agri.within + meanNL  | code + year, data = Panel_long, cluster = "code"),

  ##Elephant deaths
fepois(EDI ~ ADS2 + ADS_IN_0_1_2  + neigh + HD_dummy_0_1 +  Pop + dist + prop.Eleph250 + prop.Agri.within + meanNL | code + year, data = Panel_long, cluster = "code"))


## Getting coefficients and robust variances of Poisson regressions

Objects <- list(
  HDI = tibble(Varname =  rownames(coeftest(Regressions[[5]])), 
                  Coeff = coeftest(Regressions[[5]])[, 1],
                  Var = diag(vcov_cluster(Regressions[[5]], "code"))),
  
  EDI = tibble(Varname =  rownames(coeftest(Regressions[[6]])), 
                  Coeff = coeftest(Regressions[[6]])[, 1],
                  Var = diag(vcov_cluster(Regressions[[6]], "code"))))

##Calculating IRR = e^b and its se = sqrt(e^b*var*e^b) = IRR*se

Objects <- Map(function(x) {
    mutate(x, IRR = exp(Coeff),
           )}, Objects)

Objects <- Map(function(x) {
    mutate(x, IRR_se = IRR*sqrt(Var),
           )}, Objects)

Objects <- Map(function(x) {
    mutate(x, IRR_z = IRR/IRR_se,
           )}, Objects)


Regressions[[5]]$coefficients <- (as.vector(Objects$HDI[, 4]))[[1]] 
Regressions[[5]]$se <- (as.vector(Objects$HDI[, 5]))[[1]] 


Regressions[[6]]$coefficients <- (as.vector(Objects$EDI[, 4]))[[1]]
Regressions[[6]]$se <- (as.vector(Objects$EDI[, 5]))[[1]]

##Constructing VCOV matrices for regressions

Vcov <- list(
  vcov_conley(Regressions[[1]], cutoff = 25),
  vcov_conley(Regressions[[2]], cutoff = 25),
  vcov_conley(Regressions[[3]], cutoff = 25),
  vcov_conley(Regressions[[4]], cutoff = 25),
  (diag(exp(c(Objects$HDI$Coeff))))%*%vcov_cluster(Regressions[[5]], cluster = "code")%*%(diag(exp(c(Objects$HDI$Coeff)))),
  (diag(exp(c(Objects$EDI$Coeff))))%*%vcov_cluster(Regressions[[6]], cluster = "code")%*%(diag(exp(c(Objects$EDI$Coeff))))) 



##Table: Linear
etable(list(Regressions[[1]], Regressions[[2]]), title = "Effects of Anti-Depredation Squads on mortality due to conflict (additional controls)", vcov = list(Vcov[[1]], Vcov[[2]]), fitstat = c("n", "aic", "bic"), style.tex = style.tex(main = "base"), replace = T, signif.code = c("***"=0.01, "**"=0.05, "*"=0.10), drop.section = c("fixef"), extralines = list("Baseline" = c(Baseline_avg_012$HDI_scaled, Baseline_avg_012$EDI_scaled), "Fixed effects" = c("Two-way", "Two-way"), "Standard Errors" = c("Conley", "Conley")), dict = c("Pop" = "Population", "ADS2" = "ADS", "dist" = "dist_eleph_habitat", "neigh" = "ADS_neighbor", "HDI_scaled" = "HDI", "EDI_scaled" = "EDI", "HD_dummy_0_1" = "HD_0_1"),  label = "tbl:all_TWFE", file = "C:/Users/tanay/downloads/Appendix_all regressors_TWFEs.tex")

##Table: Poisson + YFE
etable(list(Regressions[[3]],Regressions[[4]], Regressions[[5]], Regressions[[6]]), title = "Effects of Anti-Depredation Squads on mortality due to conflict (additional controls)", vcov = list(Vcov[[3]], Vcov[[4]], Vcov[[5]], Vcov[[6]]), fitstat = c("n"), style.tex = style.tex(main = "base"), replace = T, signif.code = c("***"=0.01, "**"=0.05, "*"=0.10), drop.section = c("fixef"), extralines = list("Baseline" = c(Baseline_avg_012$HDI_scaled, Baseline_avg_012$EDI_scaled, Baseline_avg_012$HDI, Baseline_avg_012$EDI), "Fixed effects" = c("Year", "Year","Two-way","Two-way"), "Standard Errors" = c("Conley", "Conley", "Clustered", "Clustered")), dict = c("Pop" = "Population", "ADS2" = "ADS", "dist" = "dist_eleph_habitat", "neigh" = "ADS_neighbor", "HDI_scaled" = "HDI", "EDI_scaled" = "EDI", "HD_dummy_0_1" = "HD_0_1",  "ADS_nbd_ever" = "ADS_nbr_ever"), label = "tbl:all_YFE_Poi", file = "C:/Users/tanay/downloads/Appendix_all regressors_Poisson_YFE.tex")


```


#Appendix Regressions - interaction

```{r}
Regressions2 <- list(
    ##Human deaths
  
feols(HDI_scaled ~  ADS2 + ADS_IN_0_1_2 + neigh + ADS_interaction + Pop + dist + prop.Eleph250 + prop.Agri.within + meanNL | code + year, data = Panel_long, conley(25)),

   ##Elephant deaths

feols(EDI_scaled ~ ADS2 + ADS_IN_0_1_2 + neigh + ADS_interaction  +  HD_dummy_0_1 + Pop + dist + prop.Eleph250 + prop.Agri.within + meanNL | code + year, data = Panel_long, conley(25)))


##Table
etable(Regressions2, title = "Effects of ADS on mortality: Robustness checks", fitstat = c("n", "aic", "bic"), style.tex = style.tex(main = "base"), replace = T, signif.code = c("***"=0.01, "**"=0.05, "*"=0.10), drop.section = c("fixef"), extralines = list("Baseline" = c(Baseline_avg_012$HDI_scaled, Baseline_avg_012$EDI_scaled), "Fixed effects" = c("Two-way","Two-way"), "Standard Errors" = c("Conley", "Conley")), dict = c("Pop" = "Population", "ADS2" = "ADS", "dist" = "dist_eleph_habitat", "neigh" = "ADS_neighbor", "HD_dummy_0_1" = "HD_0_1", "EDI_scaled" = "ED", "HDI_scaled" = "HD"), label = "tbl:interaction", file = "C:/Users/91941/Desktop/Raghav/CECFEE/Elephant_Paper/Regression Tables/Feb10_24/Appendix tbl2_interaction.tex")


```


##Plot
```{r}

HDI_ADS <- Regressions2[[1]]$coefficients["ADS2"]
HDI_ADS_Var <- (Regressions2[[1]]$se["ADS2"])^2 
HDI_interaction <- Regressions2[[1]]$coefficients["ADS_interaction"]
HDI_interaction_Var <- (Regressions2[[1]]$se["ADS_interaction"])^2
HDI_covar <- vcov_conley(Regressions2[[1]], cutoff = 25)["ADS2", "ADS_interaction"]

EDI_ADS <- Regressions2[[2]]$coefficients["ADS2"]
EDI_ADS_Var <- (Regressions2[[2]]$se["ADS2"])^2 
EDI_interaction <- Regressions2[[2]]$coefficients["ADS_interaction"]
EDI_interaction_Var <- (Regressions2[[2]]$se["ADS_interaction"])^2
EDI_covar <- vcov_conley(Regressions2[[2]], cutoff = 25)["ADS2", "ADS_interaction"]

Interaction_table <- tibble()

Plot <- tibble(year = 2005:2016, t = year - 2005) %>% mutate(HDI_effect = HDI_ADS + t*HDI_interaction, HDI_se = sqrt(HDI_ADS_Var + t^2*HDI_interaction_Var + t*HDI_covar), EDI_effect = EDI_ADS + t*EDI_interaction, EDI_se = sqrt(EDI_ADS_Var + t^2*EDI_interaction_Var + t*EDI_covar)) %>% 
                                   mutate(HDI_CI_L = HDI_effect - 1.96*HDI_se,
                                          HDI_CI_H = HDI_effect + 1.96*HDI_se,
                                          EDI_CI_L = EDI_effect - 1.96*EDI_se,
                                          EDI_CI_H = EDI_effect + 1.96*EDI_se) 


Plot_HDI <- ggplot(Plot, aes(x = year, y = HDI_effect, ymin = HDI_CI_L, ymax = HDI_CI_H)) +
  geom_pointrange() + geom_hline(yintercept = 0) + theme_clean() +
  labs(
    x = "Year", y = "Effect on human deaths", title = "Temporal effects of ADS")


Plot_EDI <- ggplot(Plot, aes(x = year, y = EDI_effect, ymin = EDI_CI_L, ymax = EDI_CI_H)) +
  geom_pointrange() + geom_hline(yintercept = 0) + theme_clean() +
  labs(
    x = "Year", y = "Effect on elephant deaths")

ggsave(Plot_HDI, filename = "C:/Users/91941/Desktop/Raghav/CECFEE/Elephant_Paper/Plots/Temporal effects_HDI.png")

ggsave(Plot_EDI, filename = "C:/Users/91941/Desktop/Raghav/CECFEE/Elephant_Paper/Plots/Temporal effects_EDI.png")

plot_temporal <- ggpubr::ggarrange(Plot_HDI, Plot_EDI, nrow = 2, ncol = 1)

ggsave(plot_temporal, filename = "C:/Users/91941/Desktop/Raghav/CECFEE/Elephant_Paper/Plots/Temporal effects.png")

```



